%% Preliminary
clear;clc;
load('IRFs.mat');

%% Collect the results
MatSize     =   size(IRFs);
KeyResults  =   cell(size(IRFs));
for ii=1:prod(MatSize)
    [rr,cc]         =   ind2sub(MatSize,ii);
    KeyResults{ii}  =   SubFun_KeyResults(IRFs{rr,cc});
end
save('KeyResults.mat','KeyResults');

%% Data for Plots
List_Prob_ext   =   (0.00:0.1:1.00)';
Num_Prob_ext    =   length(List_Prob_ext);

List_Fin_AdjCost=   [0,0.01,0.05,0.1,0.2,0.4,0.8,2,5,10,50,100,1000];
Num_Fin_AdjCost =   length(List_Fin_AdjCost);
XAxis       =   log(1+List_Fin_AdjCost)-log(1+0.8);
YAxis       =   List_Prob_ext;

VarInfo     =   struct('C_Agg_M_ext',[1,2],'C_dom_M_ext',[2,2],'C_ext_M_ext',[3,2],...
                       'DevUIP_M_ext',[4,2],...
                       'C_Std_M_ext',[5,2],'C_GapByFinInt_M_ext',[7,2],...
                       'EffReturn_M_ext',[9,2],...
                       'C_Agg_Y_H',[1,3],'C_dom_Y_H',[2,3],'C_ext_Y_H',[3,3],...
                       'DevUIP_Y_H',[4,3],...
                       'C_Std_Y_H',[5,3],'C_GapByFinInt_Y_H',[7,3],...
                       'EffReturn_Y_H',[9,3]);
ZAxis       =   struct();
VarList     =   fieldnames(VarInfo);

for ii=1:length(VarList)
    vv          =   VarList{ii};
    Ind_Data    =   VarInfo.(vv);

    ZAxis.(vv)  =   zeros(size(KeyResults));
    for rr=1:size(KeyResults,1)
        for cc=1:size(KeyResults,2)
            ZAxis.(vv)(rr,cc)   =   KeyResults{rr,cc}(Ind_Data(1),Ind_Data(2));
        end
    end
end

%% Plots: Iso-AggC, StdC
%--------------------------------------------------------------------------
% Data for Plots
%--------------------------------------------------------------------------
Level_Baseline  =   1.6037e-04*100;
Levels      =   VecFun_MakeLogDiffGrid(Level_Baseline*0.5,max(ZAxis.C_Agg_M_ext(:)),10,0.25);% GR_TickValue('Unit',0.02).Setup(max(ZAxis.C_Agg_M_ext,0));
Levels      =   [Levels, Level_Baseline];
M           =   contour(XAxis,YAxis,ZAxis.C_Agg_M_ext,Levels);
close();
N_Levels    =   length(Levels);
Lines       =   cell(N_Levels,1);
N0          =   1;
Flag_Start  =   1;
N_Lines     =   0;
Temp_v      =   -inf;

while N0<=size(M,2)
    if M(1,N0)~=Temp_v
        Temp_v              =   M(1,N0);
        N_Lines             =   N_Lines+1;
        if Temp_v==Level_Baseline
            Idx_Baseline        =   N_Lines;
        end
        Lines{N_Lines}      =   struct('v',Temp_v,'n',M(2,N0),...
                                       'x',M(1,N0+(1:M(2,N0))),...
                                       'y',M(2,N0+(1:M(2,N0))));
        N0                  =   N0+M(2,N0)+1;
    else
        Lines{N_Lines}      =   struct('v',Temp_v,'n',Lines{N_Lines}.n+M(2,N0),...
                                       'x',[Lines{N_Lines}.x,M(1,N0+(1:M(2,N0)))],...
                                       'y',[Lines{N_Lines}.y,M(2,N0+(1:M(2,N0)))]);
        
        N0                  =   N0+M(2,N0)+1;
    end
end
Lines       =   Lines(1:N_Lines);

%--------------------------------------------------------------------------
% Figure Setup
%--------------------------------------------------------------------------
N_col       =   1;
N_row       =   1;
fig_size    =   [N_col,N_row].*[1,0.75]*1/2.5; 
fig_gap     =   [0.02,0.02];
fig_VMargin =   [0.15,0.02];
fig_HMargin =   [0.12,0.02]; 

fig_SubPlot =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                     fig_gap,fig_VMargin,fig_HMargin);
LEG         =   @(PH,LabelList)legend(PH,LabelList,'fontsize',8,'interpreter','latex',...
                                      'location','southeast');
XLabel      =   @(Label)xlabel(Label,'fontsize',8,'interpreter','latex');
YLabel      =   @(Label)ylabel(Label,'fontsize',8,'interpreter','latex');

AxisSetup   =   GR_Axis('TickFmt',struct('X','%g','Y','%.2g'),...
                        'FontSize',struct('X',8,'Y',8),...
                        'Exponent',struct('X',0,'Y',0),...
                        'Box','off');

Tick_X      =   @(X)GR_TickValue('Unit',4).Setup(X);
Tick_Y      =   @(Y)GR_TickValue('Symmetric',0).Setup(Y);
%--------------------------------------------------------------------------
% Iso-AggC Curves
%--------------------------------------------------------------------------
TempFun_LP  =   @(C)GR_Line('Style','-','Color',MyColor('Black',C),'Width',1.5);

% CList       =   VecFun_MakeLogDiffGrid(0.1,1,N_Lines,0.2)';
CList       =   0.05+Levels/max(Levels)*0.95;
% TempFun_LP  =   @(C)GR_Line('Style','-','Color',(1-C)*MyColor('Blue',0.5)+(C)*MyColor('Red',0.5),'Width',2);

Fig         =   FigureSetup(['Iso-AggC Curves'],fig_size);
ax          =   fig_SubPlot(1);
for ii=1:N_Lines
    if ii==Idx_Baseline
        PH  =   GR_Line('Style','-.','Color',MyColor('Gray'),'Width',2).Plot(Lines{ii}.x,Lines{ii}.y);
    else
        TempFun_LP(CList(ii)).Plot(Lines{ii}.x,Lines{ii}.y);
    end
%     plot(Lines{ii}.x,Lines{ii}.y);
    hold on;
end
AxisSetup.Setup(ax);
xlim([min(XAxis),log(1+2)-log(1+0.8)]);
ylim([min(YAxis),max(YAxis)]);
xticklabels(round(exp(xticks+log(1+0.8))-1,2));
XLabel('Portfolio adjustment cost');
YLabel('Share of fin. integrated hhs');
LEG(PH,'Baseline aggregate response');

print('-depsc','-r1000',Fig,['TableGraphs/IsoAggC_Curves']);
% colorbar;

%%
LEG         =   @(PH,LabelList)legend(PH,LabelList,'fontsize',8,'interpreter','latex',...
                                      'location','northeast');
%--------------------------------------------------------------------------
% C-Std
%--------------------------------------------------------------------------
[YY,XX]     =   ndgrid(YAxis(2:end),XAxis);
F           =   griddedInterpolant(YY,XX,ZAxis.C_Std_M_ext(2:end,:));
% AdjCost
Fig         =   FigureSetup(['C-Std Curves, along AdjCost'],fig_size);
ax          =   fig_SubPlot(1);
for ii=1:N_Lines
    if ii==Idx_Baseline
        PH  =   GR_Line('Style','-.','Color',MyColor('Gray'),'Width',2).Plot(Lines{ii}.x,F(Lines{ii}.y,Lines{ii}.x));
    else
        TempFun_LP(CList(ii)).Plot(Lines{ii}.x,F(Lines{ii}.y,Lines{ii}.x));
    end
    
    hold on;
end
AxisSetup.Setup(ax);
xlim([min(XAxis),log(1+2)-log(1+0.8)]);
ylim([1,9]);
xticklabels(round(exp(xticks+log(1+0.8))-1,2));

XLabel('Portfolio adjustment cost');
YLabel('Dispersion of consumption responses');
LEG(PH,'Baseline');

print('-depsc','-r1000',Fig,['TableGraphs/StdC_AdjCost_Curves']);

%%
% Share of Integrated HHs
Fig         =   FigureSetup(['C-Std Curves, along IntShare'],fig_size);
ax          =   fig_SubPlot(1);
TempLineList=   [1,Idx_Baseline,8];
TempMaxList =   [0.12,0.25,1];
PH          =   zeros(1,length(TempLineList));
for jj=1:length(TempLineList)
    ii      =   TempLineList(jj);
    TempInd =   [1>0,diff(Lines{ii}.y)>2*1e-2];
    yy      =   Lines{ii}.y(TempInd);
    xx      =   Lines{ii}.x(TempInd);
    pyy     =   linspace(0.0,TempMaxList(jj),100)';
    pxx     =   interp1(yy,xx,pyy);
    pzz     =   F(pyy,pxx);
    
    if ii==Idx_Baseline
%         PH  =   GR_Line('Style','-.','Color',MyColor('Blue'),'Width',2).Plot(Lines{ii}.y(TempInd), ...
%                                F(Lines{ii}.y(TempInd),Lines{ii}.x(TempInd)));
        PH(jj)  =   GR_Line('Style','-.','Color',MyColor('Gray'),'Width',2).Plot(pyy,pzz);
    else
        PH(jj)  =   TempFun_LP(CList(ii)).Plot(pyy,pzz);
    end
    hold on;
end
AxisSetup.Setup(ax);
xlim([min(YAxis),max(YAxis)]);
ylim([2,8]);
% xticklabels(round(exp(xticks+log(1+0.8))-1,2));
XLabel('Share of fin. integrated hhs');
YLabel('Dispersion of consumption responses');
LEG(PH,{'Low aggregate response','Baseline aggregate response','High aggregate response'});

print('-depsc','-r1000',Fig,['TableGraphs/StdC_IntShare_Curves']);
%% Plots
